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Abstract 

The new formation design strategy using simulated 
annealing (SA) optimization is presented. The SA 
algorithm is useful to survey a whole solution space of 
optimum formation, taking into account realistic 
constraints composed of contunous and discrete 
functions. It is revealed that this method is not only 
applicable for circular orbit, but also for high-elliptic 
orbit formation flying. The developed algorithm is 
first tested with a simple cart-wheel motion example, 
and then applied to the formation design for SCOPE. 
SCOPE is the next generation geomagnetotail 
observation mission planned in JAXA, utilizing a 
formation flying techonology in a high elliptic orbit. A 
distinctive and useful heuristics is found by 
investigating SA results, showing the effectiveness of 
the proposed design process. 

1. Introduction 

Japan Aerospace Exploration Agency (JAXA) is 
currently planning the next generation Earth- 
magnetosphere observation mission called “SCOPE” 
(cross-Seale Coupling in Plasma universE, Fig. I) 1 . Its 
predecessor “GEOTAIL”, launched in 1992, played a 
big role for understanding the ion behavior in the 
Earth’s magnetotail and triggered interests toward 
further microscopic scale (electrons scale) phenomena. 
Following this successful work of GEOTAIL, SCOPE 
aims at observing the Earth’s magnetotail where the 
ions and electrons interact with each other, with 5 
satellites flying in formation. To fully resolve the 
time-domain behavior and spatial distribution of the 
magnetospheric phenomena, a simultaneous 
observation by spatially distributed electro-magnetic 
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instruments is essential. 

From the formation design point of view, there are 
several requirements to be satisfied, which include not 
only the optimality of observation, but also constraints 
stemming from the operationality, the physical layout 
of the onboard components and so on. 

In the SCOPE mission, the fundamental requirements 
for formation design come from (1) observation 
optimality, (2) relative orbit determination accuracy, 
and (3) inter-satellite communication link capability 
constrained by antennae layout of each satellite. The 
(2) requirement is needed because in the SCOPE 
mission, the inter- satellite radiometric ranging and 
clock synchronization via inter- satellite 
communications are employed. (3) is needed to 
comply with the actual antennae directivity limitations 
As the first step of designing the SCOPE formation, 
one wants to obtain all the possible answers which 
satisfy (1), (2) and (3) at the same time. Here what is 
needed is surveying a whole solution space, not 
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focusing on the highly accurate local minimum, but 
rather on approximated answers which represent all 
the possible types of solutions. 

This paper proposes a formation design strategy which 
surveys a whole solution space. This method employs 
the simulated annealing (SA) algorithm to obtain the 
globally optimum solution. In this paper, the method 
is validated through analyzing the well-known 
optimum formation in low-earth orbit, and then 
applied to SCOPE formation design. Finally, based on 
the solution types obtained by SA, a heuristics to 
design the optimum SCOPE formation is introduced. 

2. Formation Optimality Index 

Before constructing the design method, this section 
introduces several indices used in the following 
formation design. 


2.1. Observation Optimality (FDOP-O) 

Consider the multi-point simultaneous observations by 
spatially distributed satellites fleet. The observation 
objective is to obtain the spatial differential coefficient 
of the observation field. The observation of /- th 
satellite, which is located at the relative position 

vector r = (x, y, z) T can be written as; 


Y=m 


(i) 


Suppose the field to be observed is linear, i.e.; 


V/C r) = 


r df df 8f V 


( 2 ) 


dx dy dz 

then from the observation of /- th and j - th satellites, the 


difference of observations 


Z tj = Y j —Y i can be 


expressed as; 

Z, l =(r J -ryvf = r ;vf ( 3 ) 

Eq.(3) implies that the quantity V/ can be calculated 
from the observation differences. To obtain 

3-dimensional V/ value, at least 3 different Z. are 


required, which means at least 4 satellites are needed 
to obtain V/ . 

Practically V/ is determined by the least-square 


method as follows; 

v /=(2>,iff'( 2>A) (4) 

where r,j is the relative position vector from /-th to j- th 
satellite. If the field is supposed to be isotropic, the 
error covariance for equation (4) can be calculated 
from the following equation; 

P = 




( 5 ) 


where 


is the error variance of the observation. 


Therefore the estimation quality can be estimated 
from the following quantity; 



( 6 ) 


which represents the dilution of precision for the 
observation quantity. (6) implies that the well-known 
tetrahedron formation is optimum, and in addition, the 
larger formation is better. By eliminating the size 
effect from (6), we finally obtain the formation 
dilution of precision for observation as follows; 


FDOPO = Jmax^. )Tr 


(LvO 


( 7 ) 


Here the index is normalized by the maximum relative 
distance among all the two satellites combinations. 


2.2. Formation Geometry Determination 
Optimality (FDOP-G) 

Consider the quality of formation geometry 
determination from the relative ranging information. 
The state vector to be estimated is the formation shape, 
which is written as; 

x_=[x l x 2 y 2 x 3 y 3 z 3 ••• y n _, z n _y 

( 8 ) 

where the Oth satellite is supposed to be at the origin, 
and the 1st satellite is supposed to be located on the 
x-axis, the 2nd satellite is on the xy-plane. Hence the 
dimension of estimation state vector becomes 3(n-2) 
for n-satellite system. 

From the relation between ranging output and relative 
position vector A/v = jr j — Y j || , the following 
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observation equation is obtained; 

dAr.. dAr.. 

SAr, = —Sr, + v - Sr , 


dr, 


dr j J 


( 9 ) 


Eq.(9) for all the satellites combinations are combined 
to form the following matrix relation; 

y = h x ( io ) 

where 

y = [£Ar 0 i SAr 02 ■■■ SAr 0n _, SAr n 


■■■ SAr ln _, ■■■ SAr n 


-in -2 3Ar n _ 2 „_ l | 


is the observation vector. The formation geometry 
estimation is now translated to minimize 

J = i(y-H x_) T (y-H x_) 


hence the dilution of precision for geometry 
determination is derived as follows; 


FDOPG = 


r~ 

/ „ X -1 


r 

1 

x' 

xl 

i 

/(«-!) 


( 12 ) 


In (12), division by n-1 is introduced to eliminate the 
satellite number dependency 
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Fig.2: Samples Distribution in E-P Plane 
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Fig.3: Correlations Between FDOPO, FDOPG and E,P 
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2.3. Formation Dynamics Determination 
Optimality (FDOP-D) 

Consider the quality of relative motion (or formation 
dynamics) determination from the relative ranging 
history information. The state vector to be estimated is 
the relative position vector (with reference to the 
inertial space), which can be written as; 

x = [jc, y, z, * 2 y 2 z 2 ••• -v i y„-i v i] T 

(13) 

Different from (8), x includes the information not only 
about the shape of the formation, but also the 
orientation of the formation with reference to the 
inertial space. Therefore the dimension of estimation 
state vector becomes 3 (n-1) for n-satellite system. 

Denoting the transition matrix of x as (t 2 , t x ) , the 

observation equation is written as y = HO r x . Then 
the formation dynamics estimation is translated to 


J = I-§(y(0- H (0®,(^ 0 )x(0)) r 

<y(0-H(0® r (f,0)x(0))A 


hence the dilution of precision for dynamics 
determination is derived as follows; 


FDOPD = 



-j) Oh, 0 ) T H( 0 r H( 0 O(/, 0 )dt 


I in- 1 ) 


(14) 

where T is the orbital period. Again division by n-1 in 
(14) is introduced to eliminate the satellite number 
dependency. 


2.4. Basic Characteristics of the Indices 

The proposed indices FDOPO, FDOPG are compared 
to the indices used in ESA’s CFUSTER mission . In 
CFUSTER, a generic evaluation of formation shape is 
done by abbreviating the formation shape to the 
pseudo-ellipsoid, and extracting the characteristic 
values “elongation” E and “planarity” P. These values 
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are derived in summary as follows. 

First, define the tensor R; 

R = Z r s r « < 3 * * * * * * * * * * * 15 ) 

The principle axes of the pseudo-ellipsoid are given 
by the eigenvectors R (n) of R. If we suppose 


R {1) > R {2) > R {3) 


(16) 


then, their square roots represent respectively the 
major, middle and minor semiaxis of the 
pseudo-ellipsoid; 

a = V/jw, b = JR™, c = JO™ (17) 


The elongation E and planarity P are defined using 
(17) as follows; 


E = \-(b/a) 
P = l-(c/b ) 


(18) 


E and P take the value between 0 and 1. The sphere 
corresponds to (E,P)=(0,0). For 4-satellite formation, 
(E,P)=(0,0) corresponds to the tetrahedron formation. 
To survey the relation between E,P and FDOPO, 
FDOPG, 1000 4-satellite formation samples are 
generated as is shown in Fig.2. Each satellite in each 
sample is generated randomly. 

The correlations between FDOPO, FDOPG and E,P 
are shown in Fig.3. Seen from these figures, the 
optimality of FDOPO and FDOPG corresponds to 
(E,P)=(0,0). It can also be said that the optimality of 
FDOPO also guarantees the FDOPG optimality. 


3. Global Optimization Using Simulated Annealing 

Algorithm 

3.1. Simulated Annealing Optimization 

To survey the whole solution space to find the 

optimum formation shapes, the simulated annealing 

(SA) optimization is employed. Some characteristics 

of SA are described below; 

• SA updates the solution by perturbing the previous 
solution randomly. 

• At first SA searches around the whole solution space, 

rather than insatiably pursuing the optimality. As the 

optimization step goes by, updates which increase 

the optimality becomes preferred more by analogy 

of a cooling process of steel “annealing”. 



Fig.4: Simulated Annealing Algorithm 


• SA is suitable for searching the discrete space, but it 
is also useful for obtaining the approximated 
optimum of the multimodal continuous space. 

• The results of SA optimization are remained to be 
refined, by investigating the proximity of the SA 
solution in detail (by local minimum search). 

The overall algorithm of SA is illustrated in Fig.4. In 
the formation optimization, the SA process becomes 
as follows; 

(1) The initial solution, which is a set of relative state 
vectors Xi=(x b y b z b v xb v yb v z i) T at a certain 
orbital phase (i.e. at apogee), is generated 
randomly. 

(2) Chose 1 satellite (/-th satellite), and regenerate 
the state vector x’i randomly. 

(3) If the objective function J(x’) is smaller than J(x), 
apply x:=x’, otherwise stochastically apply x:=x’ 
according to the process temperature T. 

(4) Terminate the algorithm when J(x) is sufficiently 
converged, otherwise update T by T:=RT and then 
go to (2). 

3.2. Validation: Cart-Wheel Motion 

To validate the SA algorithm works properly, a 
cart-wheel motion in LEO is tested. The fleet of 
satellites in a circular orbit, having relative motion 
plane 60 degree inclined from the orbital plane, 
satisfying a certain relation between relative position 
and velocity, would hold their relative geometry, while 
the whole formation is rotating one cycle per orbital 
period. It is called a “cart-wheel motion”, and it can 
be solved via Hill’s equations easily. 

If we chose the objective function so as just to keep 


4 



X X 

Fig.5: Cart-Wheel Motion Generated by SA 

the geometry of the formation, and as a result of that, 
if a cart-wheel motion is obtained by SA, it can be 
said that the SA works properly for formation 
optimization. 

The objective function is designed so as to keep the 
relative distances between all the satellites are kept to 
L for whole through the orbit; 

• / =-^§Z(h-(o Kfy (19) 

n^'2 1 i*j 


To calculate the relative orbital motion over one 
orbital period, the solution of the following relative 
orbital motion, which can be found in many papers 3 
and described briefly in appendix, is used; 


r ; +ilx i;. + 2 11 xr ; . + Qx(fixr ; .) = 


R 3 


f 

1 - 

v 


3RR 

~i R 2 


T \ 


( 20 ) 

where R is the radius, T> the orbital angular velocity, p 
the gravity constant, respectively. 

Because the relative motion must be periodic in the 
practical formation flight, the following additional 
condition must hold; 


m_ [Z ( 2+g ) 

*(0) V (l + e) 1/2 (l-e) 3/2 


( 21 ) 


y(n) = [Z (2-e) 

x(tt) ]/a 3 (l + e) 3 / 2 (l-e ) 1/2 

where the depedent variable for x,y is the true 
anomaly 0 . The result of SA optimization is shown 


in Fig.5, where the conditions are chosen as follows; 

• 3 satellites in 200km circular LEO, L=1.0 

• search space: (x,y,z)=[-l:l], (v x ,v y ,v z )=[-le3:+le3]. 

• SA setting: T 0 =100, R=0.9 

• Number of trials: 100 

In Fig.5, x is the radial direction, z the outer-plane 
direction, and y completes the right-hand system. The 
solid lines indicate the analytical solution of the 
cart-wheel motion. As the numerical solution 
distribution well aligns the analytical one, it can be 
said that the SA algorithm works properly to generate 
the optimum formation. 

4. High Elliptic Formation Optimizations 
4.1. FDOP-O Minimization Problem 

The primary objective of SCOPE mission is to 
perform multi-point simultaneous observations around 
apogee, therefore consider the following objective 
function; 

r 210° 

J = J FDOPO(6)d6 (22) 

where the reference orbit is chosen to be a high 
elliptic orbit for SCOPE((RE+3000km)x30RE). As is 
shown in section 2.4, the minimization of FDOPO 
automatically optimizes FDOPG as well. 

The result of SA optimization is shown in Fig.6, 
where the conditions are chosen as follows; 

• 4 satellites in (RE+3000km)x30RE, satO at origin. 

• Search space: (x,y,z)=[-l:l], (v x ,v y ,v z )= [-le-5: le-5]. 
Periodicity condition (21) must hold. 

• SA setting: T 0 =100, R=0.9 

• Number of trials: 100 

The coordinate system for Fig.6 is the same as that of 
Fig.5. The linear approximated analysis based on 
Eq.(20) and Eq.(Al) imply that there are 4 
geometrically equivalent formations for a given 
FDOPO and FDOPG history which are 
(x,yz)=(+,+,+), (+,+,-), (-,-,+) and (-,-,-). This sign 
selectivity has been removed and unified to (+,+,+) to 
clarify the pure formation shape dependency on 
FDOPO. Fig.6 indicates that, of course it shapes a 
tetrahedron formation, but in addition to that, there is 
a plane symmetry about yz-plane, and viewed from 
x-axis, 3 satellites (satl,2,3) forms a cocentric circle. 
In summary the FDOPO optimum formation is 
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Fig.6: FDOPO Minimum Formation Generated by SA(Left:position, Rightrvelocity distributions) 



Fig.7: Illustration of FDOPO Minimum Formation 


illustrated as Fig.7. 

4.2. SCOPE Formation Optimization 

In SCOPE mission, FDOPO is not an only thing to be 
optimized. Because of the communication 
performance limitation of the member satellite, there 
are following two constraints in real operation; 

(1) The relative distance between all two satellites 
must be within 500km to assure the inter-satellite 
communications. 

(2) Each satellite is equipped with the limited 
directivity antennae for inter-satellite 
communication. All the satellites have 
orbit-normal spin axes, and have the antenna link 
coverage of ±60deg. Therefore the relative 
location of each satellite is limited to be within 
±60deg elevation with reference to the orbital 
plane. 

Taking above (1) and (2) into account, the objective 
function for SCOPE formation optimization is written 


as follow; 


j = J FDOPO(0)dO 




A z,. 


1 / 2 
V Ax // +4>V 


-A<9 


d6 


max (Ar y (0)) 
2 min(A^.(0)) 


(23) 


where the two constraints are implemented as penalty 
functions with weighting parameters ki and k 2 . The 
second term restricts the maximum relative elevation 
between all the satellites combinations, and the third 
term tries to minimize the variation of the relative 
distance over one orbital period. 

Fig.8 is the result of SA optimization under the 
following conditions; 

• 3 satellites in 200km circular LEO, L=1.0 

• ki=1.0, k 2 =0.01, A9=40deg 

• search space: (x,y,z)=[-l:l], (v x ,v y ,v z )=[-le3:+le3]. 
Periodicity condition (21) must hold. 

• SA setting: T 0 =100, R=0.9 

• Number of trials: 100 

Fig.8 indicates that the optimum formation for 
SCOPE is not so much different in shape as that of 
Fig.6, but the velocity distribution is narrowed. It is 
understood as the effect of the two constraints, 
restricting the relative motion of the satellites in a 
whole orbital period. 


5. Heuristics for Optimum SCOPE Formation 

Based on the results obtained in section 4, the 
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Fig.8: SCOPE Optimum Formation Generated by SA (Leftrposition, Rightrvelocity distributions) 


heuristics of designing the optimum SCOPE 
formation can be constructed as follow; 

( 1 ) Periodicity condition 

All the satellites must have the same orbital 
period. 

(2) Symmetry condition 

Satl and Sat2 are on the yz-plane symmetric 
positions at apogee. Sat3 on the yz-plane at 
apogee. 

(3) Cocentric condition 

The yz-plane projection of Satl,Sat2 and Sat3 
positions are cocentric at apogee. 

If we approximate the condition (3) to the coplanar 
condition, in which SatO, Satl and Sat3 are on the 
same orbital plane, the conditions for each satellite at 
apogee can be listed as follows; 

• SatO: all 0 (on the reference orbit) 

• Satl \y lx — p , z = v z =0, y/x = l tan (a / 2) 


• Sat2: y / x = p, z = v z = 0, y / x = -/ tan (a / 2) 


• Sat3: x = x = y = z = 0, zl y — l tan (/?) 
where 


P = - 


f± (2-e) 
a 2 (l + e) 3/2 (l-e) 1/2 


and 1 is the size of the formation, a the inplane angle 
(ZSat2-0-l), (3 the outer-plane angle (elevation of 
Sat3 with reference to the orbital plane), a the 


semimajor axis. Substituting above conditions to 
Eq.(A9), we get the following formation state vector 
for each satellite at apogee; 

[Satl][Sat2] 

x — l cos(<x/2), y = ±/sin(<x/2),z = 0 


v = 


—ey 

-y(\ - e) + {\ + e) 
2-e 


ld(7r) sin {a! 2), 


v = + I0(7r)cos(a /2),v z = 0 

1 - e 


(24) 


[Sat3] 

x = 0, y = l cos /?, z = / sin [5 
v x = 0, v y = 0, v z =0 


(25) 


Originally there are 6(n-l)=18 degree-of- freedom 
design parameters for 4-satellite formation, but now 
from (24) and (25), the number of design parameters 
is reduced to only 4, that are a , 0 , y and 1. 

Because we know, from section 2, that the tetrahedron 
is the optimum shape, the unit size formation can be 
assumed as a=60deg, (3=53. 5deg, 1=1.0. Therefore 
only remaining parameter now is y . 

The simulation examples are shown in Fig. 9 and 
Fig. 10 with different y . FDOPD is 5.0 for Fig.9 
and 8.0 for Fig. 10, respectively. Seen from these 
figure, the outer-plane relative elevations are kept 
always below 60deg, and the FDOPO and FDOPG 
become small enough around apogee. It is also 
important to note that, y approximately corresponds 
to the perigee/apogee relative distance ratio. This 
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Fig.10 Heuristic Design of SCOPE Formation ( y =5.0) 

perigee/apogee relative distance ratio can be was found by investigating SA results, showing the 

approached to 1.0 at the cost of increasing FDOPO, effectiveness of the proposed design process. 

FDOPG and FDOPD. 


6. Conclusions 

The new formation design strategy using simulated 
annealing (SA) optimization was discussed. It was 
shown that SA is useful to survey a whole solution 
space of optimum formation, taking into account 
realistic constraints. It was also revealed that this 
method was not only applicable for circular orbit, but 
also for high-elliptic orbit formation flying. 

The developed algorithm was applied to SCOPE 
formation design. A distinctive and useful heuristics 
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Appendix 

The analytical solution for Eq.(20) is derived as 
follows 3 ; 

x(0) = [esin0]£/, +| 2e 2 sin 0//(0)-- 


In the same way, the state vector at apogee can be 
expressed as follows; 

x(0) = d 3 


ecos0 


(l + ecos0)~ 
j(0) =[l + ecos0]dj + [2e(l + ecos 0)//(0)]<i 2 


z(0) = 


d 2 +[-cos0]<i3 


y{Q)=[\-e\d,+ 

1 


1 - e 


[ + 1 1 Hn 0 

d 3 + 

1 

I 1 i3lll V 

Vl + ecos0 ) 

1 + ecos0_ 



*( 0 ) = - - 

l-e 

x'(0) - \_-e\d x 


(A5) 


sin0 

d 5 + 

COS0 


/(9) = - 

— + 1 

_l + ecos0_ 

_l + ecos0_ 

a 6 


l-e 


x'(0) ^[ecos©]^ + 


2e 2 cos QH(Q) + 2e 2 sin QH'(Q) + - 


esin0 2e~sin0cos0 


(l + ecos0) (l + ecos0) 

+[sin0]<i3 

/(0) =[-esmQ]d l +[2eH'(Q)-2e 2 smQH(Q) + 2e 2 cosQH\Q)~\d 2 

esinQ 


mm 


1 - e 


+ 

+ COS0 + 5- 

l + ecos0 (l + ecos0) 

m = 

cos0 esin 2 0 

rl + 

l + ecos0 (l + ecos0) 2 

a 5 + 


(l + ecos0) 
sin0 esin0cos0 


l + ecos0 (l + ecos0) 

(Al) 


where 

H(Q) = ~(\ + e 2 ) ^^--{l + e 2 )sinE + — sinEcosE 


If di is selected so that the inplane relative distance is 
preserved between apogee and perigee, it is derived 
from (A4) and (A5) that; 


[1 + 6?]^ + 
“I 

=>d x = 


1 


l-e 


1 + e 
d A 


d A = [l - e ] d x + 


1 


l-e 


(A6) 


2 U 4 


Let y denote the offset ratio of di from di 
d x = yd* 


(A7) 


(A2) 

( • )’ is 0-differentiation of ( • ), e the eccentricity and E 
the eccentric anomaly. 0-differentiation and 
t-differentiation is related by the following equation; 


• I~|7 (l + ecos6) 2 

]/a 3 (l-e 2 ) 3/2 


(A3) 


(Al) indicates that the inplane motion and outer-plane 
motion are independent of each other, and the relative 
distance history is equivalent for the sign 
combinations of (x,y,z)=(+,+,+), (+,+,-), 

The periodicity condition (21) is equivalent to d 2 =0. 
When this periodicity holds, the state vector at perigee 
can be expressed using di,. . .,d 6 as follows; 


x(0) = -d 3 
V(0) = [l + e\d x + 

z(0) = ' 1 


1 + e 


1 + e 
x'(0) = [e]d x 

1 


/(B) = 

z'(0) = 


1 + e 
1 


1 + e 


(A4) 


Substituting (A6),(A7) into (A5), we get; 
x(0) = -d 3 


y( 0)= — r — + — 
L l-e 1 + e. 

x'(0) = -eyd* 

1 


d\ 


y\ o) = 


1 + e 


-1 


and 

x(7i) = d 3 

I 

y(7r) = 


J^ + _L 

1+e l-e 
x'(7r) = —eyd x 

1 


d \ 


y\x) = - 


l-e 


-1 


d 2 


(A8) 


(A9) 


If y<l, the formation size at apogee becomes larger 
than the size at perigee, while if y>l, the size at 
apogee becomes smaller than the size at apogee. 
When y=l, the inplane distance at apogee is equal to 
that of perigee. 


9 



